Extracting falsifiable predictions from sloppy models 
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Successful predictions are among the most compelling validations of any model. Extracting falsifi- 
able predictions from nonlinear multiparameter models is complicated by the fact that such models 
are commonly sloppy, possessing sensitivities to different parameter combinations that range over 
many decades. Here we discuss how sloppiness affects the sorts of data that best constrain model 
predictions, makes linear uncertainty approximations dangerous, and introduces computational dif- 
ficulties in Monte-Carlo uncertainty analysis. We also present a useful test problem and suggest 
refinements to the standards by which models are communicated. 



Reverse engineering of biological networks entails 
working from data to models, and successful predictions 
are among the most important and compelling valida- 
tions of those models. Making useful predictions, how- 
ever, entails working back from models to data. A 
model's behavior depends both on its structure (com- 
ponents and interactions) and on its parameters (num- 
bers quantifying the structure). In biology the focus 
is on structure; parameter values themselves are gener- 
ally of little interest. To test a model structure, quan- 
titative predictions must be tempered by rigorous esti- 
mates of their uncertainties, accounting for model behav- 
ior over all sets of parameters consistent with the avail- 
able data An experimental result inconsistent with 
these uncertainties is then strong evidence that some as- 
sumption in the model structure is false. We argue that 
any valuable assessment of reverse engineering methods 
in systems biology needs to address protocols and algo- 
rithms for evaluating model prediction uncertainty. 

Here we detail how to extract falsifiable predictions 
from complex biological models. We focus on complica- 
tions introduced by sloppiness, the presence of orders of 
magnitude variation in sensitivity to different parameter 
combinations Complex biological models appear to 
be universally sloppy [2|, along with many other non- 
linear multiparameter models [3(. Sloppiness affects the 
sorts of data that best constrain model predictions 0, 
makes linear uncertainty approximations dangerous, and 
introduces computational difficulties in Monte-Carlo un- 
certainty analysis. We discuss all these issues, introduce 
a useful test system, and suggest refinements to the stan- 
dards by which models are communicated. 
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FIG. 1: Conceptual illustration and quantification of slop- 
piness. A) Contours measure the change in model behavior 
as parameters vary. The model is sensitive to stiff directions 
and insensitive to sloppy directions. Dots indicate a region 
of parameter sets consistent with available data; the region is 
naturally aligned along sloppy directions, so predictions from 
this data can have small uncertainty, even though parameter 
uncertainty is large. B) Sloppy eigenvalue spectra for models 
of: (i) growth factor signaling m , (ii) intra-receptor dynam- 
ics [3, (iii) circadian rhythms [(J, (iv) quantum Monte-Carlo 
wave functions, (v) sums of exponentials. Fitting a plane to 
data (vi), on the other hand, is not a sloppy problem. 



SLOPPINESS 



Sloppiness is illustrated conceptually in figure [I|\, 
which shows a plot in parameter space where contours 
represent surfaces of constant model behavior. The 
model is very sensitive to stiff parameter combinations 
and very insensitive to sloppy combinations. 

These sensitivities can be quantified by a cost function 



C (9) that measures the change in a model's behavior as 
parameters 9 vary. We Taylor expand C (9) about a set 
of best-fit parameters 9* yielding the Hessian matrix: 
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The eigenvectors of H are the principle axes of the el- 
lipses shown in figure [TJ^V, and the model's sensitivity to 
parameter changes along each axis is proportional to the 
square root of the corresponding eigenvalue. 

Figure |TJ3 shows eigenvalue spectra for several models 
(a subset of those in Q and 0]). Column (i) is a model of 
growth factor signaling in PC12 cells where the cost 
function measures the model deviation from a set of real 
data (48 parameters and 68 data points). Columns (ii) 
and (iii) are models of intra-receptor dynamics || and cir- 
cadian rhythms @, respectively, where the cost function 
measures the deviation from simulated data. (Here, the 
derivatives in equation [T] are taken in the logarithms of 
the biochemical parameters to reflect relative changes in 
parameter values.) In all three cases the eigenvalues span 
more than 10 6 , indicating that the models are over one 
thousand times more sensitive to some directions than 
others; that the ellipses in figure QK are over one thou- 
sand times as long as they are wide. 

Sloppiness is not restricted to biology, as illustrated by 
column (iv) in figure [1)3, which shows the eigenvalues for 
fitting parameters of a wave function for use in Quan- 
tum Monte-Carlo. Nor is sloppiness restricted to very 
complex models column (v) shows the spectrum for 
fitting the decay rates in a sum of 48 exponentials (ap- 
pendix |Aj). Not all models are sloppy, however, column 
(vi) shows the eigenvalue spectra for fitting a plane to a 
set of data (a typical form of multiple linear regression). 
The eigenvalues all have roughly the same magnitude, 
indicating that all directions in parameter space are sim- 
ilarly sensitive, so the model is not sloppy. The presence 
of sloppiness in such a diverse range of nonlinear models 
(i — v) suggests that it is a universal feature of nonlinear 
multiparameter models. 



in spite of large regions of parameter uncertainty, as long 
as that parameter uncertainty is correlated along sloppy 
directions. 

By contrast, direct parameter measurements yield un- 
corrected parameter uncertainties. For example, if 62 
were known less precisely than 61, the region of accept- 
able parameter sets would be a vertical ellipse in fig- 
ure UK- We find that generally very few bare parameter 
directions are sloppy directions, so, unless the parame- 
ters have been measured very precisely, such an ellipse 
will encompass many behaviors and predictions uncer- 
tainties will be correspondingly large. 



B. Estimation algorithms 

After optimizing to find the best-fit set of parameters, 
prediction uncertainties can be calculated by accounting 
for model behavior over the region of parameter space 
that is consistent with the data. Here we consider two 
approaches to calculate uncertainties: linearized covari- 
ance analysis (LCA) and Monte-Carlo analysis (MCA). 



1. Linear covariance analysis 

Linear covariance analysis involves two approxima- 
tions: a quadratic expansion of the cost function about 
the best-fit parameters and a linear approximation of the 
model response to parameter changes. The standard de- 
viation, a y , of the prediction y is then given by 
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dy 



(2) 



II. PREDICTION UNCERTAINTIES 

Given that sloppiness is a common feature of complex 
biological models, and nonlinear multiparameter models 
in general, we ask how sloppiness impacts making pre- 
dictions. 



LCA is computationally inexpensive, particularly in 
problems where the cost function is a sum of squared 
residuals (C (9) — J2k rk (^) 2 )- Then the Hessian matrix 
can be approximated by the Fisher Information Matrix: 
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A. Constraining parameters 

The first step in making predictions from a model is 
constraining the model's parameters. Because biological 
models are sloppy, predictions are generally much more 
efficiently constrained by collectively fitting model pa- 
rameters than by directly measuring them [2j. This can 
be understood from figure[IJ\.. Fitting parameters to data 
naturally constrains the region of parameter sets consis- 
tent with that data (indicated by the dots) to lie along 
directions to which the model is insensitive. Thus the pa- 
rameter sets encompass relatively few model behaviors 
and predictions have small uncertainties. Because the 
model is sloppy, predictions can have small uncertainties 



This approximation is very useful because for many mod- 
els (including those based on differential equations) first 
derivatives can be obtained semi-analytically, and obtain- 
ing stable finite-difference derivatives is difficult in sloppy 
problems. 



2. Monte-Carlo analysis 

Monte-Carlo analysis explicitly samples from the dis- 
tribution of parameter sets consistent with the available 
data. The posterior distribution of parameter sets given 
model M and data D is given by: 



P{6\D;M) occxp(-C(0)). 



(4) 
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is actually consistent with the data. 

The difference between LCA and MCA arises from 
nonlinearities in the cost surface and the model response. 
In particular, the sloppy directions determined by the 
Hessian matrix are exactly those directions that have 
small quadratic components in the cost. Therefore in 
these directions higher order terms are responsible for 
constraining the behavior. We observe that prior infor- 
mation on parameter ranges can mimic the effect of these 
higher order terms, allowing the LCA to produce uncer- 
tainty estimates comparable to MCA. Further, we are 
investigating the use of curvature measures of nonlincar- 
ity 0] to predict when the LCA will be inaccurate for 
sloppy models. 



FIG. 2: Comparison between linearized covariance analysis 
(LCA) and Monte-Carlo analysis. In dark grey is the one- 
standard-deviation uncertainty bound from LCA for the ac- 
tivity of Mek given NGF stimulation in our growth factor 
signaling model [j . In light grey is the corresponding Monte- 
Carlo result. The dashed line is prediction from the best-fit 
set of parameters. Note how dramatically the LCA prediction 
overestimates the uncertainty. 

After an ensemble of parameter sets has been drawn from 
this distribution, the uncertainty on any quantity can be 
estimated simply by calculating it for each member of the 
parameter ensemble and noting the resulting variation. 

Numerous algorithms exist to sample distributions 
like ((3]), but obtaining a well-converged ensemble can 
be computationally challenging. After comparing pre- 
dictions from LCA and MCA we will discuss specific al- 
gorithms and speed-ups for MCA. Importantly, however, 
once obtained for a given model and set of constraining 
data the ensemble can be used for any number of predic- 
tions. 



3. Comparison 

The approximations involved in LCA may introduce 
significant artifacts, as illustrated in figure^ which plots 
a particular prediction of our growth factor signaling 
model [ij]. LCA dramatically over-estimates the uncer- 
tainty in the prediction, particularly at early times, re- 
ducing the power of the prediction to test the model 
structure. For example, an experimental result of 3/4 
activity at 8 minutes lies outside the bounds of the well- 
converged MCA prediction and would suggest that some 
assumption in the model structure is probably incorrect. 
The LCA prediction however, would not offer this insight. 

In figure [2 LCA overestimated the uncertainty, but 
in some cases LCA may also underestimate uncertainty. 
Such underestimation may be even more damaging, as it 
may cause one to reject a model that MCA would reveal 



C. Efficient Monte-Carlo for Sloppy Models 

Although sloppiness implies that many sets of param- 
eters will be consistent with the data, those consistent 
sets remain a very small fraction of the entire parameter 
space. Directly sampling from the posterior distribution 
(equation [?]) is thus infeasible. We use Markov-Chain 
Monte-Carlo to sample the distribution, building our en- 
semble via a random walk through parameter space Q . 

For sloppy models the cost function is very stiff in some 
directions and very sloppy in others. If we take steps at 
the scale of the stiff directions, exploring the sloppy di- 
rections will be very slow, but if the step size is too large 
moves in the stiff directions will cause very few steps 
to be accepted, also slowing convergence. Thus it is vi- 
tal to use importance sampling, taking larger steps in 
some directions than others. We find that it is natural 
to scale steps using the square root of the regularized 
Fisher Information Matrix and that this gives reasonable 
acceptance probabilities. Furthermore, we find it impor- 
tant to recalculate the FIM periodically to maintain a 
reasonable chance of accepting attempted moves. This 
indicates that the cost surface is not only very narrow, 
but also substantially curved. 

Building parameter ensembles in sloppy systems can 
also challenge algorithms that solve the model equations. 
For example, the ensemble built for our differential equa- 
tion model of growth factor signaling [U explored regions 
of parameter space where the equations were very dif- 
ficult to integrate, even for well-tested stiff integrators. 
Resulting errors must be caught, and we found it useful to 
dynamically tighten integration tolerances. Significantly 
for biochemical networks, Michaelis-Menten and Hill type 
equations can become particularly difficult to numerically 
integrate when exploring parameter space because the 
"turning points" in these curves can become incredibly 
sharp and without sufficient guidance the integrator can 
overstep these points and land on a different, unphysical 
solution to the equations. 



4 



III. STANDARDS FOR THE COMPUTATIONAL 
SYSTEMS BIOLOGY COMMUNITY 

The Systems Biology Markup Language (SBML) [§] 
has emerged as a community standard for the definition 
of models of biomolecular reaction networks. By making 
our own code — SloppyCell [l(| — able to read and write 
SBML models, we have been able to apply our analyses to 
a number of published models in the literature, demon- 
strating the ubiquity of sloppiness Q. We found that 
this standard, while adequate for summarizing finished 
models, is inadequate for capturing the larger process of 
reverse engineering, in particular the derivation of model 
structure and parameter values from data. In addition, 
since biological experimentation often involves creating 
topological variants (mutants) of a wild-type, biological 
modeling requires working with multiple related networks 
which collectively define a "model". While SBML con- 
tains constructs for defining networks, it provides no data 
structures or operations for modifying or relating net- 
works. The biological reverse engineering community will 
— at some stage — need to address these deficiencies in 
current standards if community progress is to continue. 



IV. CONCLUSIONS 

Model validation is a crucial part of reverse-engineering 
biological networks, and making predictions is one of 



the most important parts of validation. To usefully test 
model structures, predictions must be accompanied by 
uncertainty estimates that account for the underlying pa- 
rameter uncertainty. Linearized covariance analysis may 
introduce important artifacts, so when possible Montc- 
Carlo analysis should be preferred, and such analysis has 
specific numerical challenges that need to be overcome. 
As the community and associated standards develop we 
anticipate that techniques to make and test predictions 
will prove increasingly important. 



APPENDIX A: FITTING EXPONENTIALS 



The fitting of a sum of exponentials is a classic prob- 
lem 11], and, as seen in figure |TJ3 it is sloppy. Fitting ex- 
ponentials has proven a very useful test problem both for 
analysis and for algorithm development. It is physically 
motivated, but computationally convenient; the cost and 
Hessian matrix are simple functions of the true decay 
rates and initial amounts. This facilitates generation of 
high quality statistics to investigate questions concern- 
ing sloppiness, such as the relationship between parame- 
ter degeneracy and the spacings between eigenvalues [1| . 
It also makes fitting exponentials a useful test case for 
comparing and refining optimization algorithms. 
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